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We calculate electrostatic potential landscapes for an external probe charge in the presence of a 
set of metallic islands. Our numerical calculation in three dimensions (3D) uses an efficient grid 
relaxation technique. The well-known relaxation algorithm for solving the Poisson equation in two 
dimensions is generalized to 3D. In addition, all charges on the system, free as well as induced 
charges, are determined accurately and self-consistently to satisfy the desired boundary conditions. 
This allows the straightforward calculation of the potential on the outer boundary using the free 
space electrostatic Green's function, as well as the calculation of the entire capacitance matrix of 
the system. Physically interesting examples of nanoscale systems are presented and analyzed. 

PACS numbers: 02.70.-c, 41.20.-q 



Introduction 

There is a need to precisely know the electrostatic 
landscape experienced by electrons in ever smaller struc- 
tures, down to the scale of scanning tunnelling micro- 
probes and single electron devices. The presence of con- 
ducting leads for manipulating and measuring local po- 
tentials influences the quantum mechanical behavior of 
electrons in a highly non-trivial manner. In polarizable 
media, for example, charged "conglomerates" which in- 
clude free as well as polarization charges in the neigh- 
borhood, behave as quasi-particles which can live for a 
comparatively long time and interact with the conducting 
leads via Coulomb interaction. Knowledge of the poten- 
tial landscape describing this interaction for the case of a 
quasi-particle, is an interesting and important element in 
the better understanding of these systems. Electrons in 
single-electron transistors 1] , or moving in the neighbor- 
hood of lithographically-defined gate arrangements [2|, 
or tunnelling through STM scanning tips, are but a few 
examples of the pervasiveness of electrostatic potentials 
in realistic structures. 

The solution of the Laplace or Poisson equation to 
obtain electrostatic potential landscapes is a well de- 
fined boundary value problem, typically requiring the 
discretization of space on a convenient grid. Relaxation 
techniques are well known 0, 0, and widely used in the 
solution of these problems, as they provide convenient 
and efficient algorithms Q. For dimensions lower than 
three, simple second order algorithms are used together 
with common speedup features such as successive overre- 
laxation (SOR) and Gauss-Seidel (GS) iteration schemes 
0, H- In three dimensions, however, due to the poor 
scaling with grid dimensions, more efficient routines are 
desirable. In this context a generalized 0(h 6 ) algorithm 
for 3D is presented in this paper, and used to calculate 
the potential landscape of several physical systems of in- 
terest, as those mentioned above. 

The boundary conditions most easily built into relax- 
ation algorithms are fixed voltage surfaces, with the volt- 
age known and provided by an external battery, for ex- 



ample. This does not apply, however, to cases with a 
floating potential, such as metallic islands which are iso- 
lated from the environment (like metallic quantum dots), 
or for open boundary problems. In the case of isolated 
islands, the value of the potential at a metallic boundary, 
even though constant, is not known. On the other hand, 
the overall charge on the island is determined at the out- 
set and can be considered to be known. The solution to 
the problem taken here is that once one has access to 
the linear relationship between the charge on the island 
and the island potential derived from the relaxation al- 
gorithm, one can invert this relationship and calculate 
the potential with every relaxation cycle such that the 
overall charge is maintained at a fixed value, e.g. zero 
for an overall neutral island. Incorporation of this idea 
in the iteration procedure yields the appropriate floating 
potentials, as we will show. 

Notice moreover that the outer boundary is open in 
general in most nano-sized geometries. If the size of the 
calculated cell could be chosen large enough, of course, 
one could assume that the potential would drop to zero 
there; however, this is operationally forbidden by the 
vast number of grid points needed in that case. The 
only feasible way is to determine the non-uniform poten- 
tial on the outer boundary self-consistently within the 
algorithm. The approach taken in this paper is that the 
knowledge of the total charge distribution (external and 
induced charges) allows for the calculation of the poten- 
tial on the outer boundary via the standard free-space 
electrostatic Green's function, l/Airr. This is a straight- 
forward if computationally expensive procedure, which 
can be speeded up remarkably by tabulating the inverse 
distances on the grid, but yields physically well-behaved 
asymptotics in all cases, as we will see. We should also 
emphasize that once the charges and potentials in the 
system are known, one can easily evaluate the capaci- 
tance matrix for the geometries of interest, regardless of 
the symmetries of the arrangement. 

The remainder of the paper describes the algorithm in 
detail in section I, while section II illustrates its use in 
several physically relevant examples. 
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FIG. 1: Visualization of averages taken on the grid: average 
(f) c as in Eq. (fTa)> is over nearest neighbors (NN); average 
(f) s as in Eq. (|lb|> is over third nearest neighbors (TNN). The 
grid points included in the sum are shown as filled symbols; 
solid lines joining these have length of two grid spacings, 2h. 



I. ALGORITHM 

Taking the standard Taylor expansion of a smooth 
function f(x,y,z) around each grid point, one defines 
the center average and the square average as follows (as- 
suming uniform grid spacing h in all three dimensions) 
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where (T)NN stands for (third) nearest neighbors and 
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with i, j, fc and r, s, t being integers. The averages in 
Eqs. are shown graphically in Fig. Notice that 
the odd-order derivatives in eqs. {Q cancel due to the 
symmetric combinations around the grid points included 
in the averages. Notice also that second (or next) near- 
est neighbors are considered in the "checkered lattice" 
sweeps of the points making the simple cubic grid. (The 
relaxation sweeps are done sequentially over the face- 
centered cubic array of neighbors which form effectively 
a dual lattice.) 

Taking the linear combination of the averages (|la|) and 
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then with a = a^D = f the overall average becomes 
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Since we are solving for the Poisson equation 

V 2 / = -£, (5) 



Eq. 



can be rewritten as 
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where (g) c is the center average for the source term at a 
given grid point. Equation JHJ), together with (|3J), serves 
as the basis for the iterative scheme: the potential / 
is relaxed to its minimum considering external sources 
g and the potential of the (third and) nearest neighbor 
sites through the linear 3D average. In order to calculate 
charges on the grid, Eq. Q is employed again and gives 
a straightforward way for calculating real (i.e. free) and 
induced charges on the grid on an equal footing. Thus 
using Poisson's equation (0, the charge contribution of 
each grid point can be calculated as follows 



qi = h 3 gi = h (-/i 2 V 2 /i) 



= yA(/i-((/))3 D +0(A 4 )). (7) 

For convenience, throughout the paper we adopt the 
following convention on units: The charge of an electron 
is taken to be 1 and the Coulomb energy is written in 
units of eV, as E = q\ [e] • <?2 [e] / (47rr); this straightfor- 
wardly implies a unit for the distance of [r] = 18.1 nm. 
In short, the units chosen are 

[energy) = eV; [charge] = e; [distance] = 18.1nm (8) 

Therefore taking these units, brings one naturally into 
the realm of nanostructures. 



A. Successive over- relaxation and iteration scheme 

The general method for successive over-relaxation 
(SOR) is described for 2D in H|| for a TV x TV array, 
and generalized here to 3D as given by 



where 
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Nj is the grid size in the jth direction, and f^ new ^ i s cal- 
culated according to Eq. JSJ. The SOR parameter uo is 
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in the range 1 < uo < 2 as required for the algorithm 
to converge. The basic idea behind uo is that if one is 
heading in the right direction (e.g. towards the solu- 
tion), why not go a bit further. An uo too large (uo > 2), 
however, results in instability of the algorithm and the re- 
laxation process overshoots and diverges. Equation (|9bj) 
was tested for different N x , N y and N z values and it was 
indeed this specific combination that gave the optimal 
value for uo (note that the window for a good uo is quite 
narrow in general). The specific structure of (|9b[) can 
be intuitively understood as follows: the SOR algorithm 
introduces perturbations in the system that propagate 
during the relaxation cycles and eventually die out if 
the grid is large enough; however, for finite grid sizes, 
the perturbations are reflected at the boundaries and so 
they interfere and pile up. In this sense, the minimum 
extension within the three grid dimensions constrains the 
optimal magnitude of uo, consistent with Eq. (|9b|). 

For further optimization of the algorithm, the Gauss- 
Seidel iteration scheme was adopted, as well as the alter- 
nating relaxation on the two checkerboard like sub-grids 
that in sum span the whole grid y, l| • The inverse dis- 
tances between the grid points were mapped into a table, 
such that the calculation of the potential in the grid is 
sped up remarkably. According to Coalson 6] , multigrid 
methods can be applied to account for the slowly con- 
verging long wavelength portions of the solution. This 
was not done here, since the variation of the potential on 
the isolated islands and especially the calculation of the 
outer boundary already introduced longer range correla- 
tions over the grid that presumably made the algorithm 
converge faster in our case. 

A note about efficiency: As we use a successive over- 
relaxation method to iterate the potential on the grid, 
the total relaxation time for this in 2D is proportional to 
~ n 3 / 2 , where n is the total number of grid points : 4], and 
is thus clearly comparable to algorithms like conjugate 
gradient. Notice also that SOR has still known improve- 
ments that may also be implemented and would thus 
make this algorithm superior to the former 4]. Our re- 
laxation over the bipartite lattices composing the simple 
cubic 3D grid preserves the spirit of the 2D algorithms, 
but obtains an accuracy of 0(/i 6 ), as discussed above. 



B. Open outer boundary 

Equation Q gives a consistent higher order recipe for 
calculating the overall charge distribution, including in- 
duced as well as the external (free) charges 0, Q , given 
via the source term g in the Poisson equation (0. Start- 
ing with an (arbitrary) initial constant potential on the 
outer boundary (OB), these values are updated every 
time the interior of the grid has relaxed down to a certain 
accuracy level e, 
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FIG. 2: Calculated interaction potential of a point charge 
(q = 1 [e]) with a finite size conducting plane. The grid con- 
figuration is shown to scale in the inset: the total block size 
is 32x32x132, indicated by the gray box lines; positions of 
the point charge are indicated by the black dots in the in- 
set, a set of locations vertically away from the metallic plane 
which is shown at the bottom by the black line (note that 
closest point-charge to island distance is 3/i, three grid spac- 
ings). Shown are two configurations: neutral isolated island 
(filled circles) and grounded island (empty circles) . Gray lines 
indicate different power law dependencies oc z a . 

and employs the tabulated inverse distance values for the 
grid points. 



II. DISCUSSION 

In the following, several applications of the algorithm 
are presented. We start with a simple test example, and 
follow with the analysis of more complex geometries. 



A. Example: point charge near a conducting island 

As an instructive example and as a test case for the al- 
gorithm, a 32 x 32 x 132 grid was setup with one square 
metallic island in the lower region, while an external 
charge is placed at different positions away from the is- 
land surface and directly over its center (see inset of Fig. 
EJ). The total dimensions of the grid in real space were 
taken to be L x = 1 (in units as per|HJ), and accordingly 
L y = 1 and L z = 4.2 for equal grid spacing. 

We calculate the interaction potential experienced by 
the point charge in the presence of the island. The results 
are shown in Fig. [21 For the case of a neutral island, close 
to the surface, the potential approaches that produced 
by the image charge of an infinite plane (~ while 
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FIG. 3: Interacting point charge (q — +1 ) with an array 
of four islands in a plane underneath the point charge: a 
64 x 64 x 16 (box) cell with open boundaries was setup; its 
dimensions are L x = L y = 1, L z = 16h = 0.25 (units as per 
®). The separation of the external charge from the island 
plane is 5h = 5/64 = 0.078. 



further away from the island, the potential approaches 
asymptotically the form of an induced quadrupole inter- 
action potential (~ 1/r 6 ) |9|- 

The calculation was done also for a grounded island. 
In this case, one obtains the limit of an infinite plane at 
short distances (~ 1/r), while further away from the is- 
land the dependence weakens to ~ 1/r 2 as expected [HI • 
As a reference, the interaction potential in the presence 
of an infinite plane is shown (solid line in Fig. EJ). The 
finite size of the island clearly reduces the interaction at 
large distances. Overall, it is interesting to see that one 
obtains exactly what is expected, but more importantly, 
the algorithm already gives useful results for rather small 
grid dimensions. 



B. Example: point charge interacting with an 
array of islands 

Lithographically, an array of conducting islands can be 
separated from a 2D electron gas by an insulating layer, 
as in the experiments with spatially modulated two- 
dimensional electron gases in semiconductors Q. Con- 
sidering the electron g as as a Fermi liquid, the interac- 
tion of a single electron (or single quasi particle) with the 
conducting islands nearby is an important element of the 
physics of this problem since this interaction will clearly 
modify the dynamical behavior of the system |l|. Fig- 
ure [3| shows the sample geometry used in this example. 

With the test charge hovering over the center of 
one island and considering the conducting islands either 
grounded or isolated (uncharged), the calculated charge 
distribution and the potential in the plane of the islands 
are shown in Fig. 21 the grounded case can be understood 
as that of an island where hopping onto and off the is- 
land is allowed via a tunnelling channel from an external 
charge reservoir. A few points should be stressed: first, 
note that the variations of the potential in the plane of 
islands in Fig. 2] for the case of isolated islands is about 
a factor 10 larger compared to the grounded case since 




FIG. 4: Potential landscape and charge distribution within 
the plane of islands for the geometry in Fig. the left (right) 
two graphs are for the case of neutral (grounded) islands, 
respectively. Notice the quite different scale in the two po- 
tentials shown in A and B. The induced charge is shown as 
surface charge in the plots C and D where for better visualiza- 
tion the top (positive values) and bottom surface (negative) 
charges of the islands are shown together in a single plot. 



the island potential is not fixed. Second, the isolated is- 
lands are indeed neutral within numerical accuracy (the 
sum over all charges Q tot is zero, within the accuracy 
provided by the convergence parameter e (see lower two 
plots in Fig. 21 total induced charge in the plane of the 
islands = all positive surface charge + all negative sur- 
face charge ~ 0.442 - 0.442 - 10" 7 - e). Third, the 
induced negative charge in the case of neutral islands is 
clearly smaller compared to the grounded case, which is 
intuitively clear, since in the neutral case the negative 
charge needs to be compensated by an equal amount of 
opposite charge, and this charge separation in the island 
costs energy. Therefore, the interaction with the exter- 
nal charge can be expected to be weaker for the isolated 
(neutral) islands, as is the case (see later, Fig.[5J). Notice 
also that in the neutral case the island corners exhibit 
accumulation of induced charges, as one would expect 
from the sharpness of the island corners. It is also inter- 
esting to observe that the induced charge never exceeds 
the external charge (in absolute value): at the maximum, 
it is just equal and opposite in sign (as for the case of the 
infinite plane). 

As the external charge was displaced horizontally at a 
certain separation above the plane of the island, one effec- 
tively scans the potential landscapes for different geome- 
tries and different boundary conditions, as shown in Fig. 
The isolated (neutral) islands show a clearly weaker 
modulation of the interaction potential. Furthermore, 
decreasing the island separation such that the gap be- 
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FIG. 5: Potential landscape of an external charge interacting 
with the array of islands. The geometry is the same as in Fig. 
El except for the separation between islands. Only the region 
of the four overlapping island corners is shown. The left two 
graphs, A and C, are for the case of isolated islands; the right 
two, B and D, for grounded islands; the upper two graphs are 
for an island separation s of lOh (= 10 grid spacings), while 
the lower two graphs are for an island separation of 3h. 

tween them is nearly closed, reduces most of the central 
modulation and smooths the potentials, as one would ex- 
pect, and can be seen by comparing the upper two graphs 
in Fig. [5] with the lower two. 

C. Mutually induced charge on STM tip 

As a final example, the induced charges where calcu- 
lated for the typical geometry of a metallic STM tip near 
a metallic array structure with the islands kept neutral; 
the geometry can be seen in panel D of Fig. The 
parabolically shaped STM tip is maintained at a fixed 
potential (V = 1) by an external source. The potential 
below the islands is weakened and shielded by the islands 
(panel A), while on the islands the potential is constant 
and is defined there by the constraint of neutrality (panel 



B). Above the islands, in the region of the tip, the po- 
tential landscape (xy-slice) has a circular plateau at the 
position of the tip with a potential of V = 1 (the po- 
tential of the tip) and smoothly decays away from the 
tip (panel C). Panel D indicates the charge distribution 
on the islands and on the tip; the charge on the tip is 
positive overall: Q% v = +4.01 is the charge on the por- 
tion of the tip shown in the absence of the islands, and 
it gains a bit more of charge (+0.181) through the in- 
teraction with the metallic islands. The island with the 
tip right on top of it (right rear in inset), is the most 
polarized of all. The corner in the center is negatively 
charged, while the transition through the white region on 
the island surfaces in panel D indicates that the surface 
charge switches sign there, and that the outer region of 
the islands are positively charged. In addition, the lower 
island surfaces are also mostly positive, as one would ex- 
pect. This compensates the negative charge induced by 
the tip and guarantees the neutrality of the islands (see 
caption for explicit numbers). 



III. SUMMARY 

In summary, the electrostatic potential of complex 
metallic arrangements were calculated on a three dimen- 
sional grid with an 0(h 6 ) algorithm. The algorithm pre- 
sented here is a generalization of the relaxation tech- 
niques common in 2D systems, properly set up to provide 
accurate and efficient calculations. The approach allows 
the study of arbitrary geometries and boundary condi- 
tions, as well as the self-consistent calculation of free and 
induced charges. This information, in turn, allows the 
calculation of the capacitance matrix of the system. Sev- 
eral examples illustrate the reliability and usefulness of 
the algorithm for obtaining potential landscapes of inter- 
est. 
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FIG. 6: Potential and charge distribution for a conducting 
tip at potential V = 1 above four isolated (neutral) islands; 
the distance between the plane of islands and the tip is 5/i. 
Panels A to C show x?/-slices in the potential distribution. 
Panel A: potential just under the plane of islands (z = —Ah 
with respect to the plane of the islands); panel B: potential 
right on the plane of islands (z = Oh); and panel C: potential 
on a plane above the tip (z = 25h). Panel D is a contour 
plot of the charge distribution on the tip and on the overall 
neutral islands (Qi s l t arids = +0.281 + (-0.281) ~ 10~ 10 , i.e. 
nicely converged). The extra charge on the tip due to the 
presence of the islands is +0.181. 
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in [l(j. 
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ductor can be obtained by minimizing the energy of the 
probe charge in the presence of the charged conductor, 
E = Q 2 /2C + qQ/z, where q is the probe charge, Q than 
on the conductor of capacitance C, and z is their sepa- 
ration. The minimal energy yields Q — —qC/z, so that 
the potential on q is then qQ / z — —q 2 C/z 2 . 



